%This file generates a case of a 20x20 Cov-matrix with up to 10 sub blocks.
%10 sampled sequences of same length are generated for use in the
%expected value objective function. Noise is also added to each block.
%The performance is measured for different number of blocks and different
%sample lengths.

clear, close all
NumOfStocks = 20;
NumOfSamples=300;
MeanSample=10;

%[CovBlocks, Combinations, SampleSizes, Sigma, StockReturns]=GetCovBlocks(0);
figure(1)
load data_10_blocks.mat

%Create a 10 Sampled sequences of numbers. Maximum sample number = 5000.
Samples=randn(NumOfStocks,NumOfSamples);
Samples=Sigma^0.5*Samples;

%Create noise
SampleNoise=randn(NumOfStocks,NumOfSamples,MeanSample);

error=[];
errorNoise=[];

for iterBlocks=2:5
NumOfBlocks=2*iterBlocks;

%Define the covariance blocks for this iteration instance
CovBlocksIter=CovBlocks(1:NumOfBlocks);
CombinationsIter=Combinations(1:NumOfBlocks);


%Define number of samples per block
SampleLengths=NumOfSamples*ones(1,NumOfBlocks);

%Define weights according to sample length, make them sum 1.
w=zeros(1,NumOfBlocks);
for i=1:NumOfBlocks
    w(1,i)=SampleLengths(1,i)/sum(SampleLengths);
end

    for k=1:NumOfBlocks
        CovBlocksSampled{k}=cov(Samples(CombinationsIter{1,k},1:SampleLengths(1,k))');
    end

    cvx_begin
        cvx_quiet(true); 
        variable Sigma_hat(NumOfStocks, NumOfStocks) symmetric;
    
        %Define objective function
        f = 0;
        for a=1:10
            for q=1:NumOfBlocks
                f = f + w(1,q)*norm(Sigma_hat(CombinationsIter{q}, CombinationsIter{q})-CovBlocksSampled{q},'fro');           
            end
        end
        minimize (f)
    
        subject to
        
        Sigma_hat == semidefinite(NumOfStocks)
     cvx_end
    
%Now test the result using frobenius norm

error(iterBlocks-1)=norm(Sigma-Sigma_hat,'fro');
index=1;

    NoiseMagnitudes=0.3*[0 2 6 10 15];
   %Loop over different noise levels
    for iterNoise=1:length(NoiseMagnitudes)
        NoiseMag=NoiseMagnitudes(iterNoise);
        %Disturb the samples with noise
        for i=1:MeanSample
            SamplesIter(:,:,i)=Samples+NoiseMag*eye(NumOfStocks,NumOfStocks)*SampleNoise(:,:,i);
        end
        
        %Create the noisy Covariance blocks
        CovBlocksSampledNoisy = cell(MeanSample,NumOfBlocks);
        for j=1:MeanSample
        for k=1:NumOfBlocks
            CovBlocksSampledNoisy{j,k}=cov(SamplesIter(CombinationsIter{1,k},1:SampleLengths(1,k),j)');
        end
        end

            cvx_begin
                cvx_quiet(true); 
                variable Sigma_hatWithNoise(NumOfStocks, NumOfStocks) symmetric;

                %Define objective function
                f = 0;
                for a=1:MeanSample
                    for q=1:NumOfBlocks
                        f = f + w(1,q)*norm(Sigma_hatWithNoise(CombinationsIter{q}, CombinationsIter{q})-CovBlocksSampledNoisy{a,q},'fro');           
                    end
                end
                minimize (f)

                subject to

                Sigma_hatWithNoise == semidefinite(NumOfStocks)
             cvx_end

        %Now test the result using frobenius norm
        errorNoise(index,iterBlocks-1)=norm(Sigma-Sigma_hatWithNoise,'fro');
        index=index+1;
        iterNoise
        norm(Sigma-Sigma_hatWithNoise,'fro')
    end

    AbsoluteError(:,iterBlocks-1)=errorNoise(:,iterBlocks-1) - error(iterBlocks-1);
end
    
figure(1)
    hold on
    
    plot([0:10],errorNoise(:,1),':k')
    plot([0:10],errorNoise(:,2),'-*k')
    plot([0:10],errorNoise(:,3),'--k')
    plot([0:10],errorNoise(:,4),'-.k')

    legend('4 Blocks', '6 Blocks', '8 Blocks', '10 Blocks')
    title('Estimation error')
    xlabel('Level of noise on submatrices')
    ylabel('norm(S-Shat)')

figure(2)
    hold on
    
    plot([0:10],AbsoluteError(:,1),':k')
    plot([0:10],AbsoluteError(:,2),'--ok')
    plot([0:10],AbsoluteError(:,3),'--k')
    plot([0:10],AbsoluteError(:,4),'-.k')

    legend('4 Blocks', '6 Blocks', '8 Blocks', '10 Blocks')
    title('Increase in estimation error due to noise')
    xlabel('Level of noise on submatrices')
    ylabel('increase in error')
